// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Copyright (C) April 27, 2010 MARIN - All rights reserved - http://www.marin.nl
// -----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File     : IterCudaSolver.h
// -----------------------------------------------------------------------------

#ifndef ITERCUDASOLVER_H
#define	ITERCUDASOLVER_H

#include "Solver.h"
#include "Array2D.h"

#define ITERCUDASOLVER_BASE Solver

// C++ flaw: declarations of private variable types needed
namespace iterative_cuda
{
    template <typename, typename> class cpu_sparse_csr_matrix;
    template <typename, typename> class gpu_sparse_pkt_matrix;
    template <typename, typename> class gpu_vector;
    template <typename>           class diagonal_preconditioner;
}

// -----------------------------------------------------------------------------
// class IterCudaSolver
// -----------------------------------------------------------------------------
/// iterative CUDA solver by Andreas Klöckner
/** http://mathema.tician.de/software/iterative-cuda */
// -----------------------------------------------------------------------------
class IterCudaSolver : public ITERCUDASOLVER_BASE
{
    typedef ITERCUDASOLVER_BASE base;

public:

    IterCudaSolver(Array2D<REAL> *pcc,
                   Array2D<REAL> *pcs,
                   Array2D<REAL> *pcw,
                   int Nx1,
                   int Nx2,
                   REAL tolerance);

    virtual ~IterCudaSolver();

    virtual void solve(Array2D<REAL> *pB,
                       Array2D<REAL> *pX);

private:

    /// @cond nevercall
    IterCudaSolver(const IterCudaSolver& m);
    IterCudaSolver & operator =(const IterCudaSolver& m);
    /// @endcond

    void preconditioner(Array2D<REAL>& matc,
                        Array2D<REAL>& mats,
                        Array2D<REAL>& matw);

    void makeMatrixMarketFile(Array2D<REAL>& matc,
                              Array2D<REAL>& mats,
                              Array2D<REAL>& matw,
                              const char* file);

    typedef iterative_cuda::cpu_sparse_csr_matrix<REAL, int> CpuMatrixType;
    typedef iterative_cuda::gpu_sparse_pkt_matrix<REAL, int> GpuMatrixType;
    typedef iterative_cuda::gpu_vector<REAL, int>            GpuVectorType;
    typedef iterative_cuda::diagonal_preconditioner<GpuVectorType> Preconditioner;

    CpuMatrixType*  m_pCpuMatrix;      ///< matrix on CPU
    GpuMatrixType*  m_pGpuMatrix;      ///< matrix on GPU
    GpuVectorType*  m_pInvDiagGpu;     ///< inverted diagonal (?)
    GpuVectorType*  m_pPermInvDiagGpu; ///< permuted inverted diagonal (?)
    Preconditioner* m_pDiagPre;        ///< preconditioner (?)
    REAL*           m_pB;              ///< right hand vector
    REAL*           m_pX;              ///< result vector
    unsigned int    m_n;               ///< row count of sparse matrix
};

inline Solver* NewIterCudaSolver(Array2D<REAL> *pmatC, ///< central matrix
                                 Array2D<REAL> *pmatS, ///< south matrix
                                 Array2D<REAL> *pmatW, ///< west matrix
                                 int N, ///< width
                                 int M, ///< height
                                 REAL tolerance) ///< solver tolerance
{
    return new IterCudaSolver(pmatC, pmatS, pmatW, N, M, tolerance);
}

#endif // ITERCUDASOLVER_H

